1. Field of the Invention
This invention relates to a diffusion simulation method by a computer simulation of a process of manufacture of a semiconductor device, and more particularly to a diffusion simulation method of diffusion of an impurity into the inside of a device.
2. Description of the Related Art
In a computer simulation of impurity thermal diffusion, which is one of steps of manufacture of a semiconductor device, a domain to be analyzed is divided into a mesh and a diffusion equation is discretized at each mesh point, whereafter the diffusion equation is linearly converted into simultaneous linear equations by a Newton's method to determine a solution to the diffusion equation as disclosed in Ryo Dan, "Process Device Simulation Technique", pp.26-28.
In a condition of a low temperature and a high concentration in which an impurity is a solid solution, the diffusion of the impurity cannot be described in a simple diffusion equation which only includes a total concentration of the impurity as a variable, but requires such a diffusion model which takes clustering into consideration as is described, for example, by expression (12) on page 398 and expressions (19) to (21) on page 399 of S. M. Sze, "VLSI TECHNOLOGY", 1983. If the diffusion model is interpreted more generally, then it is described by the following three equations. ##EQU1##
Here, expression (1) is a diffusion equation which includes a diffusion term by a movable impurity and an electric field drift term by an electrically active impurity, and expressions (2) and (3) are rate equations for defining the concentration of the movable impurity and the concentration of the electrically active impurity, respectively. C is the total concentration of the impurities, C.sub.m the movable impurity concentration, and C.sub.a the electrically active impurity concentration. D is the diffusion constant of the impurity, and .mu. is the drift mobility of the impurity by the electric field. Z is the valency of the electrically active impurity having the positive or negative sign added thereto, and .psi. is the electrostatic potential by which an external driving electric field is formed. Further, f(C, C.sub.m) is an expression which represents a rate at which the impurity can be rendered movable, and g(C, C.sub.a) is an expression which represents a rate in which the impurity is electrically activated. While it is considered that C.sub.m =C.sub.a stands for almost all impurities, also an impurity whose diffusion is suitably modeled as C=C.sub.m &gt;&gt;C.sub.a like diffusion of phosphor doped in a high concentration in silicon is available. Therefore, it is more general to handle f(C, C.sub.m) of expression (2) and g(C, C.sub.a) of expression (3) as separate expressions from each other.
When it is tried to actually solve expressions (1) to (3) on a computer, it is required to perform discretization in the space and time domains to convert them into algebraic equations which can be solved by the computer. For the space discretization, a control volume method is used wherein a mesh is produced in an object analysis domain, and the Gauss' theorem: ##EQU2## is applied around each mesh point to remove diverging operators, whereafter slope operators are replaced by finite differences. In equation (4) above, .OMEGA. represents the domain assigned to each mesh point, and .differential..OMEGA. represents an outer surface portion which forms the boundary of the region assigned to each mesh point. Further, A represents an arbitrary aimed for vector amount, and n represents a normal vector to the outer surface. More particularly, for example, to such a two-dimensional mesh structure as shown in FIG. 1, expression (1) is applied after it is transformed as ##EQU3## signifies that a sum is calculated between aimed for mesh point i and an adjacent mesh point coupled to mesh point i by mesh branch ij, and S.sub.ij represents that portion of domain .OMEGA..sub.i assigned to the aimed for mesh point which is a contribution by mesh branch ij. .omega..sub.ij represents that portion of outer periphery .differential..OMEGA..sub.i of .OMEGA..sub.i which is a contribution of mesh branch ij. C.sub.i represents the total impurity concentration at mesh point i. J.sub.ij represents an impurity flux which flows out from mesh point i to mesh point j, and D.sub.ij and .mu..sub.ij represent the diffusion constant and the electric field mobility at the middle point on mesh branch ij, respectively. l.sub.ij represents the length of mesh branch ii, and C.sub.mi and .psi..sub.i represent the movable impurity concentration and the electrostatic potential at mesh point i, respectively. C.sub.au represents the movable impurity concentration at a mesh point which corresponds to the upstream side of the drift flow by the electric field, and a countermeasure for an upstream difference as just mentioned is taken in order to assure numerical stability.
Meanwhile, the discretization in the time domain is performed by replacement of the time differentiation term by a difference by using a finite time step size. For evaluation points of time of the terms other than the time differentiation term, an implicit solution such as a backward Euler's method which uses an evaluation value at the present point of time of analysis is employed frequently for stability in numerical calculation. In the present example, the following equations obtained by applying a backward Euler's method to expressions (5), (6), (2) and (3) discretized already in terms of the space so as to be discretized also in terms of the time are used as equations to be solved finally: ##EQU4## where t.sub.n represents the present point of time, and t.sub.n-1 represents the directly preceding point of time of analysis. .DELTA.t.sub.n represents time step size t.sub.n -t.sub.n-1 at present. Where the equations given above are non-linear with regard to C.sub.i (t.sub.n), C.sub.mi (t.sub.n), C.sub.ni (t.sub.n) and so forth, they are usually solved after a Newton's method is applied to convert them into linear equations.
The conventional solution of a diffusion equation described above is illustrated as a flow chart in FIG. 2.
First, in step 301 of FIG. 2, a mesh is produced in an internal domain of a semiconductor device whose impurity diffusion is to be simulated. Then in step 302, initial concentration values of individual diffusion species are set on each of the mesh points produced in step 301. In step 303, a point of time of analysis is initialized, and in step 304, the analysis time point is incremented by one. In step 305, an impurity diffusion equation, a movable impurity concentration definition equation and an electrically active concentration definition equation are set as simultaneous equations and solved collectively to determine an impurity concentration. Finally in step 306, it is checked whether or not an analysis end time point set in advance is reached. If the analysis end time point is reached, then the simulation is ended, but if the analysis end time point is not reached, then the control returns to step 304.
It is a problem of the technique proposed by the conventional example that a diffusion equation must be set as simultaneous equations and solved together with another equation for determining a movable impurity concentration and a further equation for determining an electrically active impurity concentration. When simultaneous linear equations converted into linear equations by a Newton's method are solved by a direct method, the order of calculation cost is 0(M.sup.2) where the total number of equations is M, and where the number of mesh points is N, since M=3N, the calculation amount is 0(9N.sup.2).
Consequently, when compared with a calculation amount of 3.times.0(N.sup.2) required to solve the three equations individually, a calculation time as long as approximately three times is required. Also with regard to a storage area, approximately 3 times is required to store the three equations simultaneously.